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Abstract. - We study the asymptotic behaviour of the Bak, Tang, Wiesenfeld sandpile 
automata as a closed system with fixed energy. We explore the full range of energies char- 
acterizing the active phase. The model exhibits strong non-ergodic features by settling into 
limit-cycles whose period depends on the energy and initial conditions. The asymptotic activity 
p a (topplings density) shows, as a function of energy density £, a devil's staircase behaviour 
defining a symmetric energy interval-set over which also the period lengths remain constant. 
The properties of C,-p a phase diagram can be traced back to the basic symmetries underlying 
the model's dynamics. 



Sandpile automata and the associated self-organized critical scenario have been widely 
used to study the occurrence of avalanche behavior in nature £Q. These automata are, in 
general, driven-dissipative systems in which matter (sand) or energy is externally added to 
the system and dissipated by the dynamics itself. Eventually, the input and output balance 
produces a stationary state with highly fluctuating bursts of activity, that in the limit of an 
infinitely slow input (time scale separation) is self-similar 4 5 6] . This condition is customarily 
embedded in the original sandpile model introduced by Bak, Tang and Wiesenfeld (BTW) 0, 
and the Manna model PI, which define a deterministic and stochastic relaxation dynamics, 
respectively. 

Recently, fixed energy sandpiles (FES) automata, which share the same microscopic dy- 
namics of the corresponding BTW sandpiles without external driving and dissipation, have 
been studied UHSJOHEI- The energy density or sand is therefore a conserved quantity that 
acts as a tuning parameter. FES present an absorbing state phase transition (APT) jTT at a 
value of the energy density which is identical to the stationary energy density of the driven- 
dissipative version of the model jH] ■ FES with stochastic dynamics (Manna dynamics) belong 
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to the universality class of APT with a conserved field [121113) and the FES with deterministic 
BTW dynamics defines a different critical behavior with anomalies associated to non-ergodic 
effects 8 , similarly detected in the original driven model |14U15| . 

Interestingly, FES automata allows the study of the full phase diagram of sandpile au- 
tomata. In particular, the overcritical phase is relevant with respect to experimental situations 
since many avalanche phenomena are naturally poised in this regime. Indeed, relations among 
sandpile automata, charge density waves and linear interface models have been established 
and prompt to interesting behavior for the higher total energy properties |16lll7lHHlll9| . 

In this Letter we study the FES in order to provide a more detailed characterization of 
the model over the whole energy range. The density of active sites shows a step-like behavior 
for increasing energies, with constant plateaus in correspondence of energy intervals which 
form a hierarchical and symmetric interval set. The model shows strong non-ergodic features, 
and after a transient relaxes onto periodic orbits which depend upon energy and the initial 
conditions. Both the period lengths and transient times to reach the periodic ortbits remain 
constant onto the same energy intervals charactering the plateaus of the activity behavior. We 
tested the robustness of the observed behavior by looking at the scaling of activity and period 
lengths with the system size. As a preliminary understanding of some of the features observed 
in numerical simulations, we discuss analytically some basic symmetry properties of the BTW 
toppling dynamics that account for the symmetry observed in the activity behavior and the 
corresponding energy intervals. The present results might also be relevant to intermittency 
and predictability issues in slowly driven BTW automata [20112 1| and could provide new 
insights for the study of toppling invariants and recurrent states [22) . 

Here we consider the original BTW model in the square lattice with periodic boundary 
conditions. The configuration is specified by giving the energy, z^, at each site. The energy may 
take integer values, and is nonnegative in all cases. Each active site, i.e., with (integer) energy 
greater than or equal to the activity threshold zth ( z i > z th — 4), topples and redistributes 
its energy following the updating rules z% — > z% — zth, and Zj — > Zj + 1 at each of the 4 
nearest neighbors of i. The BTW dynamics with parallel updating (all active sites topple at 
each update) is completely deterministic and Abelian; i.e. the order in which active sites are 
updated is irrelevant in the generation of the final (inactive) configuration [21] . 

In the FES, the energy density £ = L~ 2 zi is a conserved quantity since no external 
driving is present and periodic conditions are assumed at the boundaries of the lattice (no 
dissipation) . The value of £ is fixed by the initial condition which is generated by distributing 
£L 2 particles randomly among the L 2 sites. Once particles have been placed, the dynamics 
begins. If after some time the system falls into a configuration with no active sites, the 
dynamics is permanently frozen, i.e., the system has reached an absorbing configuration. By 
varying £, FES show a phase transition separating an absorbing phase (in which the system 
always encounters an absorbing configuration), from an active phase with sustained activity. 
This has been assumed to be a continuous phase transition, at which the system shows critical 
behavior. The order parameter is the stationary average density of active sites p a , which 
equals zero for £ < £ c . Previous studies focused on the behavior close to the critical point in 
order to characterize the scaling behavior p a ~ (£ — Cc)' 3 , for C > Cc E! Numerical results, 
however, pointed out possible failures of simple scaling hypothesis and nonergodic behavior as 
also observed in slowly driven simulations |14II15| . These features contradict to the standard 
scaling observed in stochastic sandpile models [Sj. 

In order to exploit the effect of the deterministic BTW dynamics we consider the active 
sites density behavior for the whole range of £ > Cc', i-e. the active phase. Since for any finite 
size, the number of possible configurations is finite and the dynamics is deterministic, after a 
transient the system enters a periodic orbit, visiting recursively a finite set of configurations. 



F. BAGNOLI et al: SHORT PERIOD ATTRACTORS IN FIXED ENERGY SANDPILES 



3 




2.0 3.0 4.0 5.0 

c 

Fig. 1 - Stationary density p a = nr/T vs. energy density for lattice sizes: L — 50 (square), L = 100 
(open circles), L = 200 (triangles). Data are averaged over N = 200 random initial conditions. The 
solid line curve represents the activity p a for the stochastic Manna model (zth = 4) with size L — 100. 
Inset: enlargement of the region 4.23 < C < 4.26 (box) showing that a smaller plateaus are observed 
on finer scales. 



In such a case, p a can be computed as the ratio of the total number of topplings at a site in 
a period and the length of the period itself. In Fig^we plot the activity density p a in the 
stationary state as a function of the energy (. 

The activity has a stairway structure with constant activity plateaus over energy intervals 
distributed symmetrically with respect to the energy value £ = 3.5. The step- like behavior is 
in contrast with the smooth and regular curve obtained in the case of a stochastic models in 
which, for each toppling event, energy grains are redistributed at random between neighbors 
|3] . The value of p a corresponding to each value of £ is the average over at least N — 200 initial 
conditions (IC). In the activity plateaus we do not observe a dependence of the stationary 
activity with respect to IC. For energies in the central and larger plateau, for instance, all 
initial conditions evolve to an orbit with period 2 and a single toppling per site per period for 
large lattice size. Correspondingly, the observed value of activity is delta-like distributed at 
p a = 1/2. On the contrary, for energies outside of the plateaus, we find that different IC's may 
generate different values of the stationary p a which exhibits a moderate dispersion around the 
mean reported in the plot of Fig. ^ The mean value is not recovered in individual runs, 
signalling the non-ergodic properties already noticed in Ref. [Hj- It is worth remarking that 
this scenario presents striking analogies with complex phase-locking plateaus characterized by 
devil's staircases found in numerous non- linear driven systems |2*3*] . 

A confirmation of the non-ergodic nature of the BTW model is provided by measuring, as 
a function of energy, the length of the periodic cycle where the same sequence of topplings is 
executed over and over (periodic orbit). The period lengths depend on the initial condition 
and on the system's energy. In Fig. El we plot the average period length (T) as a function of 
£ for a system of size L = 100. In the same figure p a (suitably rescaled) is shown to better 
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Fig. 2 - Behavior of average period (T) as a function of energy density for system size L — 100. The 
average is performed over at least 10 3 randomly chosen initial conditions. Dashed line indicates the 
corresponding stationary p a after a suitable rescaling for a direct comparison. 



appreciate the very same correspondence of the plateau structures. On the plateaus the cycles 
period assumes extremely low values, such as T = 2,4,6. As for the activity p ai the period 
lengths are delta distributed in correspondence of the plateaus; i.e. all IC define the same 
T . Outside the energy plateaus, different IC generate different orbits defining a nontrivial 
period distribution. It is interesting to note that a similar irregular behavior of periods is 
found in other deterministic dissipative systems with continuous state variables which exhibit 
unpredictable transient dynamics but admit asymptotic periodic states as attractors |25j . The 
robustness of the present picture has been tested through simulations of systems of increasing 
size. In FigQwe report data for system sizes L = 50, 100 and 200. The discrete nature of the 
system allows to discriminate energies on a finer scale (1/L 2 ), and the structure of the energy 
intervals defined by activity plateaus reveals progressively a reacher structure [21] • With data 
at hand we cannot say whether, in the infinite size limit, the system is characterized by an 
energy range made exclusively by infinitesimal intervals and correspondingly, the curve p a {C) 
is discontinuous. In this case also the transition inactive-active would be discontinuous, ruling 
out the possibility even to define the critical exponent j3 characterizing the transition. 

It is also worth discussing the non-ergodic effect with respect to IC. This feature does 
not disappear for increasing sizes, even though a single specific value for the period and the 
number of topplings per period (and therefore for p a ) emerges with overwhelming frequency. 
As an illustration of the basic phenomenology, we report in Fig. 01 the normalized histograms 
of periods T collected from a sample of 10 3 random IC, at different sizes and C — 4.25. The 
results indicate that, as the size L grows, each distribution seems to becomes more and more 
peaked around a given period value Tm {Tm — 8 for Q — 4.25) and the inset of Fig.[3]shows the 
increase of the histogram maximum Pm (peak) with L. This could be consistent with a period 
selection scenario in which, in the thermodynamic limit, all IC select a single period length. 
Finite size effects, however, are still too large to conclude definitely that such a mechanism 
occurs for L — > oo. 
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Fig. 3 - Histograms of periods collected in the stationary states at the energy density £ = 4.25 and 
different system sizes: L — 20 (circles), L — 50 (dashed) ,L = fOO (dot-dashed) and L = 200 (solid). 
The maxima P m of the histogram, corresponding to T — 8, lie outside of the graph. The inset shows 
the increasing of P m with the system size L, this numerical evidence seems to support the view that 
period selection becomes stronger as the size increases. 

A first understanding of the behavior of the fixed energy BTW model can be achieved 
analytically by exploiting the symmetry contained in the deterministic dynamical rules. The 
non-ergodicity with respect to the IC can be traced back to conserved invariants of the dynam- 
ics. For instance, the initial configuration in which all grains are concentrated in the central 
site of the lattice will necessarily evolve preserving its central symmetry and asymmetric 
configurations will never be visited. Other quantities conserved by the dynamics can be con- 
structed mimicking the center of mass and the angular momentum of an isolated Hamiltonian 
system. A straightforward calculation shows that 

MX = J2 Xty xz x<y 

MY = J2 x , y yz*,y (1) 

MXY = Y, x , y {x 2 -y 2 )z x , y 

are dynamical invariants if taken modulus the linear size of the lattice, with x and y indicating 
the integer coordinates of the lattice sites with respect to the origin (x = 0, y = 0). The phase 
space relative to a given value of the energy is therefore partitioned according to the values 
assumed by conserved quantities |27j . A detailed and comprehensive discussion on toppling 
invariants can be found in Ref. |2b|. It is worth noticing, that the previous considerations hold 
for any finite system but it is not possible to generalize them to the thermodynamic limit. 

A striking symmetry property of the model concerns the distribution of toppling per site 
within any given periodic orbit. In order to exploit this property let rij be the number of times 
that the site i topples during a period. Since after a period the initial configuration is restored, 
each site has given to its neighbors as many grains as it has received, which translates in: 
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Am=J2 n h ( 2 ) 

where j € nn denote the sum over the nearest neighbors of the site i. Let us suppose that 
not all rii are the same and define i* the site with the minimum number of topplings in a 
period . The quantity pi = rij — ri{* satisfies the same equation as the n's and p^ = is 
the minimum p. Since pi* is the sum of four p's, that, by definition are nonnegative we have 
that all p's in the neighbors of i* are zero. Iterating the argument, the entire lattice can be 
covered, and one obtains that all p's are zero and therefore all n's are the same. This results 
immediately implies that any site topples exactly the same number tit of times over a given 
period of length T and eventually leads to the interesting result 

Pa = "jT- ( 3 ) 

Finally, we can focus on the evident symmetries of the activity behavior p a {0 and the corre- 
sponding energy intervals. The overall symmetric structure of the energy intervals centered 
at £ = 3.5 can be accounted for by the following argument. Let us consider two energy fields 
Zi{t) and Zj'(t), that at time t are related as Zi(t) — 7 — z'^t), but otherwise arbitrary. It can 
be easily verified that the above relation holds at all subsequent times. In fact, the evolution 
equation for z^s can be written as 

Zi (t + 1) = Zi (t) + M*) > 4] - 40[zi(t) > 4], (4) 

and by substituting Z{(t) — ► 7 — z[{t) we obtain 

z t (t + 1) = 7 - z[(t) + e[zj (t) < 3] - AQ[zi(t) < 3]. (5) 

By using the relation Q[g < 3] = 1 — 9[<? > 4] we readily obtain 

«i(t+l) =7-z' i (t + l), (6) 

This show recursively that if at a given time t the two fields are related by Zi(t) = 7 — z[(t) 
they are similarly related at all subsequent times. In addition if a site topples in one of the 
two lattices, the corresponding site in the other won't, and vice versa leading to p a {zi) = 
1 — Pa(z'i)- The full account of the observed symmetry under transformation £ — > 7 — £ and 
Pa(C) ~ * 1 ~ Pa{7 — C) would require that an equal proportion of IC settle into configurations 
related by z[ — 7 — Zi. While it is not possible to show analytically the latter statement, the 
symmetry is recovered exactly in the diagram of Fig. ^ 

Despite these simple symmetry arguments account for several features of the BTW model, 
many other issues remain unsettled. In particular, the activity and period length plateaus 
extending over the energy intervals do not find a simple analytical explanations. As well, 
the continuous or discontinuous nature of the p a (C) behavior cannot be discriminated by the 
simple analytical arguments provided here. 

In summary, we presented a numerical and analytical study of the BTW automaton with 
fixed energy. We find strong non-ergodic features related to the deterministic dynamics of the 
model in the whole energy range. Some features of the model can be understood in terms of 
the basic symmetries of the dynamics. A full rationalization of the present results could help 
to understand the scaling anomalies and the universality class of deterministic self-organized 
critical sandpile models. 
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